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The simulation of low-temperature properties 
of many-body systems remains one of the ma- 
jor challenges in theoretical [IHH] and experi- 
mental [4— 6j quantum information science. We 
present, and demonstrate experimentally, a uni- 
versal cooling method which is applicable to any 
physical system that can be simulated by a quan- 
tum computer [7149]. This method allows us to 
distill and eliminate hot components of quantum 
states, i.e., a quantum Maxwell's demon [10] . The 
experimental implementation is realized with a 
quantum-optical network, and the results are in 
full agreement with theoretical predictions (with 
fidelity higher than 0.978). These results open 
a new path for simulating low-temperature prop- 
erties of physical and chemical systems that are 
intractable with classical methods. 

From quantum field theories pQ to high-Tc supercon- 
ductivity [2] and chemical reactions [3], quantum simu- 
lation [7]-[9] provides a unique opportunity for both the- 
orists and experimentalists to explore a domain of sci- 
ence that goes beyond the applicability of any known 
classical computing method. Modern low-temperature 
physics has advanced primarily due to the development 
of efficient cooling methods [TTJ [12]. The same is also 
true for quantum information science. Physical cooling 
can drive quantum states towards the physical ground 
states, which are pure states. However, in the con- 
text of quantum simulation, pure states are not neces- 
sarily "cooler" than mixed states. For example, pure 
states with an equal superposition of all eigenstates cor- 
respond to infinite-temperature states [13j [14] . In order 
to achieve cooling for quantum simulation, in general, 
one not only needs to employ physical cooling to avoid 
decoherence [15], but also be able to prepare states that 
have low entropy in the basis of the eigenstates of the 
Hamiltonian of the system being simulated. 

An attractive approach for achieving cooling for spins 
(or qubits) is the heat-bath algorithmic cooling (HBAC) 
method [16] . The main idea of HBAC is to reduce the 
entropy of spin systems by unevenly distributing more 
entropy to one of the spins that can release the excess 
entropy to a heat bath through thermalization. The fea- 
sibility of HBAC has been demonstrated experimentally 
with NMR technology [4 . However, HBAC is not a uni- 
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FIG. 1: (Color online). Overview of the cooling method, (a) 
Logic diagram of the feedback cooling system. The cooling 
module produces two outcomes, correlated with heating and 
cooling. The measurement result can be mapped into the po- 
sition x of a ID random walker; when the walker goes beyond 
its starting position x = to the negative position x = — 1, 
the particle is either rejected or recycled, (b) The quantum 
circuit diagram of the cooling module. It includes a controlled 
evolution for time t and an energy bias parameter 7 for opti- 
mizing the performance. 



versal method for cooling quantum many-body systems; 
it is primarily employed for preparing polarized spins 
as initial states, or ancilla qubits, for quantum compu- 
tation. Furthermore, the temperature of the environ- 
ment imposes a fundamental limit for the cooling perfor- 
mance [T7] . 

A more general approach p~8] [19] for quantum cool- 
ing is to engineer a dissipative open-system dynamics to 
drive quantum states to the ground states of the sim- 
ulated systems. This dissipative open-system (DOS) 
approach relaxes two criteria in the standard circuit 
model [20] for quantum information processing, namely 
specific initial-state preparation and full-unitary gate de- 
composition. As the underlying dynamics is based on 
Lindblad master equations, this approach is shown to be 
robust [21] against simulation errors. A proof-of-principle 
demonstration has been experimentally performed using 
trapped ions [5 . Despite the advantages the DOS ap- 
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FIG. 2: (Color online). Basic principle of the cooling method, 
(a) The relative change in the population of the output state 
post-conditioned with a "cooling" measurement result de- 
pends on the eigenenergy Ek : the lower the energy, the higher 
the gain. Details are summarized in (b). (see also Eq. [I]and 
method summary). 



proach can provide, the range of application of it is lim- 
ited; practically, it is restricted to a class of system in- 
volving the so-called frustration- free Hamiltonians [19], 
where the ground state should simultaneously minimize 
all local terms in the Hamiltonian. This condition is not 
satisfied by the majority of physical Hamiltonians. 

To overcome the challenges, in this work, we present a 
universal cooling method applicable to cooling any phys- 
ical system that can be simulated by a quantum com- 
puter. Our method shares some features in the HBAC 
and the DOS approach, but does not suffer from the lim- 
itations of them. Essentially, we replace the requirement 
of the existence of a heat bath in HBAC by an ensemble 
of quantum states; thermalization is replaced by swap- 
ping with one of the states in the ensemble. As evident 
by our numerical analysis (see [22 J, our method does not 
suffer from similar restrictions in HBAC for cooling spin- 
1/2 particles [17 . Furthermore, we adopted a feedback 
approach to result in an effective quantum operation for 
achieving cooling for quantum simulation. In DOS, each 
of the quantum jump operators encodes only partial in- 
formation of the Hamiltonian to be simulated JTSJ [19]; 
in our case, the jump operators contain full information 
about the Hamiltonian. This property makes our cooling 
method universally applicable to any physical Hamilto- 
nian. 

To be more specific, for any given Hamiltonian H s and 
any input state, pure state \ipi n ) or mixed state p^ n , the 



method guarantees that the energy E = Tr (H s p out ) of 
output state p out is less than (or equal to) that of the in- 
put state. Importantly, compared with DOS, our method 
does not increase the qubit resource required for cool- 
ing; other than the system qubits, it requires only one 
extra ancilla qubit, which is the same as in the DOS 
method. This not only significantly reduces computa- 
tional resources but also makes the method more ro- 
bust against noise and errors [21 . Therefore, our cool- 
ing method is as realizable as DOS with currently avail- 
able technologies; we experimentally performed a proof- 
of-principle demonstration with quantum-optics, which 
is reported below. 

The core component of the method is the 'cooling 
module' formed by a simple quantum circuit shown 
in Fig. [l] It consists of an ancilla qubit initialized 
into the |0) state, a pair of Hadamard gates H = 
(l/x/2) [(|0> + |1» <0| + (|0) - |1» <1|], a local phase gate 
7 = |0) (0| — ie 11 |1) (1|, and a controlled unitary opera- 
tion U = exp (— iH s t), where H s is the Hamiltonian to 
be simulated. We note that the unitary time evolution 
for most of the physical Hamiltonians can be efficiently 
simulated with a quantum computer [7]. For any input 
state l^in), the quantum circuit produces the following 
state [221: 



A_|^ n )|0)+A + |^ n )|l) 



(1) 



where A± = (I =b ie ll U)/2. Next, the ancilla qubit is 
measured in the computational basis {|0) , |1)}. Similar 
to the working mechanism of DOS [TBI 033 EQ, the 
non-unitary operators A± play similar roles as quantum 
jump operators for describing the state evolution. More- 
over, we found that (see [22]) the energy of the resulting 
state A_ \ipi n ) (A_|_ \i/Ji n )), after nomalization of the vec- 
tor norms, is lower (higher) than that of the input state 
\ipin). Consequently, the cooling module probabilistically 
projects any input state |^ n ) (or p in for mixed states), 
into either a higher-energy state A + \ipi n ) (ancilla output 
|1)) or a lower-energy state A_ \ipi n ) (ancilla output |0)), 
with respect to the Hamiltonian H s being simulated (see 
Fig, [g^i and [2}}, and [22 for a step-by-step derivation). 
Eq. [^represents the key element in our cooling method. 

In the next step, the resulting state A± \ipi n ) will ei- 
ther be sent to another cooling module for further cool- 
ing, or recycled/rejected (see method), depending on the 
measurement outcome. We formulated this procedure 
as a feedback-control loop as depicted in Fig. [T^t. As 
governed by the laws of quantum mechanics, the mea- 
surement outcomes of the ancilla qubit are intrinsically 
random. We model the cooling process as a ID random 
walk by keeping track of the sequence of the measure- 
ment outcomes. A net cooling effect can be achieved 
by implementing the feedback system that employs the 
following mapping: consider a random walker starts at 
position x = 0. If the measurement outcome is |0), then 
x is increased by 1, and decreased by 1 for |1). When 
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FIG. 3: (Color online). Experimental details, (a) The input photon state is prepared by a polarization beam splitter (PBS), a 
half- wave plate (HWP) and a quarter- wave plate (QWP). The cooling module is a polarization-dependent Sagnac interferometer. 
Split by the PBS, the horizontal and vertical components of a polarized photon propagate in opposite paths within the 
interferometer. Two HWPs operate on the corresponding polarization components. The quartz plate compensator (PC) is 
used to compensate for the relative phase in the interferometer. These elements simulate the system Hamiltonian. Both paths 
are then recombined on the same PBS. The photons at output ports and 1 represent the cases of cooling and heating, 
respectively. The photon proceeds to the next cooling module depending on the value of the position x of the random walker. 
The final quantum state is reconstructed using quantum state tomography, with the measurement bases defined by a QWP, 
HWP and PBS. The photon is detected by a single photon detector (SPD) equipped with a 3 nm interference filter (IF). When 
the simulated Hamiltonian is a x , a HWP with an angle of 22.5° and a tiltable QWP with an horizontal angle implement the 
a x transformation at the input and output ports, (b) Sketch of the random walker evolution. 



the random walker goes beyond the origin to the nega- 
tive side (x < 0), it is either discarded ("evaporative" 
strategy) or recycled ( "recycling" strategy) where the 
system is reset to the the initial state |^ n ) (see [22] for 
more details). In implementing the evaporative strat- 
egy, we assumed an initial ensemble of identical states 
(red particles in Fig. [T^i). When the random walker po- 
sition is set to x = — 1, then the system is rejected, and 
we restart the experiment with another member in the 
ensemble. Just like ordinary evaporative cooling, this re- 
sults in a decrease of overall energy but at a cost of loss 
of particles (blue particles in Fig. [T^i). For the recycling 
strategy, it is applicable to a single system. We assume 
the same initial state can be prepared efficiently. When 
the random walker position is set to x = — 1, we simply 
reset the quantum state of the system to the same initial 
state. For x / 1, both strategies work in exactly the 
same way. 

We note that the recycling strategy resembles the way 
of releasing entropy through thermalization in HBAC. 



The main difference is that we correlate high-entropy 
states with the position (x < 0) of the random walker 
instead of including a pre-assigned qubit as heat sink in 
HBAC. The performance and trade-off of the two strate- 
gies will be compared and explained in the experimental 
section. In both cases, the random walker undergoes a 
diffusion toward the x > direction (see Fig. J3Jd). A 
detailed scaling analysis with Monte Carlo simulation is 
provided in [22] . 

In fact, the cooling module is effectively a realization 
of the thought experiment of Maxwell's demon [10]: an 
intelligent being who is able to separate hot and cool 
particles, which decreases the system's entropy without 
any work being done on it. In our setup, ancilla qubit s 
are constantly refreshed, and they therefore serve as an 
entropy sink. Each measurement outcome of the cool- 
ing module reveals partial information about the eigen- 
states of H s . If we repeatedly apply the cooling mod- 
ule to the same system, the process is asymptotically 
close to a quantum non-demolition measurement [23] in 
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the eigenstate basis of H s (see [22 for details). In other 
words, entropy (in the energy basis) is removed from the 
system by measurements, instead of employing an ex- 
ternal heat bath as in the heat-bath algorithmic cooling 
method [4j [16] . 

Below, we report a proof-of-principle demonstration of 
the quantum cooling method with an all-optical setup 
(see e.g. [24 for a recent review of photonic quantum 
simulation). We built a network of optical elements to 
connect four cooling modules (see Fig. |3^i) to implement 
a multiple-step quantum cooling. This setup is specifi- 
cally designed for cooling the polarization degrees of free- 
dom of a single photon (see [22] for requirements for a 
scalable implementation). The simulated quantum sys- 
tem has only two levels, i.e., a qubit. To illustrate the 
flexibility of the setup, two different Hamiltonians H s 
were simulated, namely Pauli's matrices a z = ( o -i ) 
and a x = ( i J ) • We emphasize that although the Hamil- 
tonians to be simulated are relatively simple, the ground 
states are not assumed to be known. The other specifica- 
tions of the simulation are as follows: the logical Hilbert 
space of the system was encoded in the photonic polariza- 
tion {\H} , \ V)} from a single photon source, where \H) 
and \V) respectively refer to the horizontally and verti- 
cally polarized state of a photon. The which-path degree 
of freedom {|0) , |1)} was employed as the ancilla qubit in 
the quantum circuit (see Eq. [I]) . The evolution operator 
U (t) = exp (—iH s t) was simulated for a fixed period of 
t = 7r/2. An energy offset, defined as 



7 + 7T/2 , 



(2) 



where 7 is defined in Fig. [TJd, was used as an adjustable 
parameter for characterizing the cooling efficiency. The 
goal of the experiment is to minimize the mean energy 
(a z ) (or (cr x )) for any given initial state. 

Each of the cooling modules were realized by a 
polarization-dependent Sagnac interferometer (circled by 
dotted lines in Fig. |3|i); similar, but not identical, struc- 
tures [25j [26] have been used for tasks related to quan- 
tum information processing. The quantum state \ipi n ) 
of the incident photon was prepared by a series of opti- 
cal elements shown in Fig. [3^i. The optical components 
inside the cooling module result in quantum logical op- 
erations [27] which are equivalent to that described in 
Eq. [l] (see [22] for details) . The photon leaving from one 
of the two output ports of the cooling module is in a su- 
perposition state (see Eq. [I]) containing both outcomes 
of heating |1) (red arrow) and cooling |0) (blue arrow). 
Subsequent detection of the two paths corresponds to a 
projective measurement in the ancilla qubit space, and 
will result in one of the two output states A± |^ n ). 

For the first cooling module, a single step of cooling is 
achieved if the measurement outcome is |0). To further 
extract energy from the system, multiple steps of cool- 
ing events are needed. To proceed, four identical cooling 
modules are connected by optical fibers. The network 
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FIG. 4: (Color online). Experimental results, (a) Mean en- 
ergies for 3 cooling steps of the evaporative (non-recycling) 
strategy for different 0. Different maps correspond to theory 
and experiment for Hamiltonians a z (Z) and a x (X), and dif- 
ferent initial states (step 0). (b) Theoretical and experimental 
mean energy and ground state probability after the 3rd cool- 
ing step as a function of the initial ground-state overlap |/3| 2 . 
The evaporative strategy simulated the Hamiltonian a z with 
energy bias angle = 10°. (c) Theoretical and experimen- 
tal mean energy and ground state probability after the 3rd 
cooling step as a function of 0. The initial input state was 
(2/V5)|e) + (1/V5)| ff ) (M 2 :|/3| 2 =4:l). 



of cooling modules (as shown in Fig. |3^l) allows us to 
implement both the evaporative and recycling strategies. 
The implementation was guided by following the random- 
walker diagram shown in Fig. [3}d. Specifically, for any 
event where the photon emerges from the cooling (heat- 
ing) exit of a cooling module, the random walk position 
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FIG. 5: (Color online). Experimental results. In (a-h) 
we compare the tradeoff between the evaporative and re- 
cycling cooling strategies for Hamiltonian a z , initial state 
(|e) + \g))/V%, and bias angle = 10°. (a) Theory and ex- 
periment for the number n of copies obtained out of N initial 
copies at consecutive steps. The total probability of the recy- 
cling algorithm is set to 1, correcting for experimental loses, 
(b) Experimental results for the corresponding mean energies, 
(c-e) Theoretical results for the mean energies for the evap- 
orative strategy, the recycling strategy, and their difference. 
Different bars correspond to different energy bias angles and 
cooling steps, (f-h) Corresponding experimental predictions, 
(i) and (j) are the corresponding experimental and revised fi- 
delities. Error bars correspond to the statistical distribution 
of experimental measurements. 



x is increased (decreased) by 1. Across the boundary 
x = — 1, for a range of parameters, the walkers have en- 
ergy (H s ) higher than that of the initial state. As in 
evaporative cooling, separating these hot walkers from 
the system necessarily lowers the average energy of the 
system. 

The experimental results with multiple steps of the 
cooling simulation are summarized in Fig. [4j The color- 
coded energy maps of Fig. [4^i correspond to the evap- 
orative strategy. To characterize the efficiency of the 
cooling method, we considered different initial input 
states a\e) + f)\g) for Hamiltonians H s = cr z and cr x , and 



different proportions \a\ 2 : \f3\ 2 between the populations 
of the ground \g) and excited states \e). We remark that 
(i) the cooling method does not require the knowledge of 
the ground states, and (ii) the theoretical values of the 
mean energy are the same for both Hamiltonians. The 
first row in each color coded energy map (step 0) rep- 
resents the mean energy of the initial state. The other 
rows in an energy map correspond to the mean energy 
after each step of cooling. The columns of the energy 
maps correspond to different choices of the energy bias 
= 7 + 7r/2 (cf. Eq. [2| in the cooling module. For an 
angle = 0°, the mean energy remains constant. For 
other angle settings (0 ^ 0°), the mean energy decreases 
with increasing cooling steps. Note that in the evapora- 
tive strategy, after the first step, the output states corre- 
sponding to heating events (x — —1) are discarded. This 
is reflected in the fact that for all angles ^ 0° the mean 
energy decreases (output port of the first cooling mod- 
ule). However, at the second step, the random walkers 
jump to either x = and x = 2 and no walker is rejected 
(output ports 00 and 01, see also Fig. ^jp). Consequently, 
the mean energies after the first and second steps are the 
same (see also Fig. [5^). 

We systematically probed the mean energy (H s ) and 
the ground-state probability P g = (g\ p exp \g) of the out- 
put states corresponding to the cooling events after three 
steps, Figs. [4]3 and [4^. The mean energy decreases with 
increasing ground-state probability |/3| 2 and energy bias 
0. The output state approaches the ground state when 
is close to 90°, in good agreement with the theoreti- 
cal prediction. In our experiment, the statistics of each 
count are considered to follow a Poisson distribution and 
the error bars are estimated from the standard deviations 
of the values calculated by the Monte Carlo method. We 
compared the trade-off between the strategies with and 
without recycling of "hot" copies. Figures^ a-h) show 
the theoretical and experimental results with initial in- 
put state (|e) + \g))/\/2, Hamiltonian a z , and energy bias 
angle = 10°. It is clear that the mean energy of the 
evaporative strategy decreases faster than that of the re- 
cycling strategy at the cost of a lower total yield. Fi- 
nally, the experimental fidelities are always larger than 
0.983 and the fidelities after post-processing (see meth- 
ods summary) are even better (compare black vs red data 
points in Fig. [5} and [5]). 

To conclude, our simulation results suggest that 
cooling for quantum simulation can be experimentally 
achieved with high fidelity. Due to the modest experi- 
mental resource requirement, even the fully scalable ver- 
sion (see [22]) of this method can potentially be imple- 
mented with many of the currently available quantum 
computing architectures, such as trapped ions, supercon- 
ducting devices and NMR systems. Furthermore, as the 
solution to classical [28] or quantum [29] computational 
problems can be encoded in the ground state of certain 
Hamiltonian, our cooling method can potentially have 
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applications beyond quantum simulation. 



METHODS 

Here we justify the claim in Eq. [I] that for any quan- 
tum state \ipi n ) that is not an eigenstate of H S1 the 
state A_ (A+ |V>in))> after normalization, has en- 

ergy E — (H 3 ) lower (higher) than that of \ipi n ). To 
proceed, we performed a eigenvector expansion, which 
gives \ip in ) = ^Z k c k |e/c), where |efc)'s are the eigenvec- 
tors of H a , i.e., H s \e k ) = E k \e k ), and c k = (e k \ip in ) . 
Consider A± \ip in ) = ^ k c k (l ± ie~ l(t)k ) |e fc ), where </> k = 
E k t — 7. Note that the square norm of the amplitudes 
are \c k \ 2 (1 ± sm(/) k ). It means that apart from an over- 
all constant, each of the weight \c k \ 2 of the eigenvectors 
are scaled by a factor (1 ± sin </>&). In Fig. within the 
range \(j) k \ < 7r/2, the function 1 — sin^ is monotonically 
decreasing. In other words, the lower the eigen-energy is, 
the higher the weight gains, which results in a lower av- 
erage energy, i.e., cooling. A similar argument applies for 
the case of heating. 

Here we explain the difference between the "evapora- 
tive" strategy and the "recycling" strategy. In imple- 
menting the evaporative strategy, we assumed an initial 
ensemble of identical states (red particles in Fig. [lji) . 
When the random walker position is set to x = — 1, then 
the system is rejected, and we restart the experiment 
with another member in the ensemble. Just like ordinary 
evaporative cooling, this results in a decrease of overall 
energy but at a cost of loss of particles (blue particles in 
Fig. |TJl) . For the recycling strategy, it is applicable to a 
single system. We assume the same initial state can be 
prepared efficiently. When the random walker position 
is set to x = — 1, we simply reset the quantum state of 
the system to the same initial state. For x ^ 1, both 
strategies work in exactly the same way. 

The revised fidelities in Fig. [4]i and [4^ are obtained by 
mapping each experimental density matrix (obtained us- 
ing maximum likelihood) to its eigenvector with highest 
eigenvalue. This corrects isotropic (depolarizing) noise. 
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1: BASIC PRINCIPLE OF THE SIMULATED COOLING ALGORITHM 



Here we explain in detail the basic idea of the cooling algorithm implemented in this experiment. As shown in 
Fig. [TJd, the core of the quantum circuit (cooling module) consists of four components: (a) two Hadamard gates 

H = ^(|0> + |1»<0| + ^(|0)-|1»(1| , (3) 

applied to the ancilla qubit in the beginning and at the end of the quantum circuit, (b) a local phase gate, 

R z ( 7 ) = |0><0|-^|1><1| , (4) 

where the parameter 7 plays a role in determining the overall efficiency of the cooling performance of the algorithm, 
(c) a two-qubit controlled unitary operation, 

l(g)|0)(0| + [/(g)|l)(l| , (5) 

where 1 is the identity operator, and 

U(t) =e~ iHst (6) 

is the time-evolution operator for the system. Here H s is the Hamiltonian of the system, and the energy is defined as 
the quantum expectation value of H s . The parameter t also determines the overall efficiency of the algorithm. 

For any given initial state (it works equally well for mixed states), |^ n ), the algorithm starts with the product 
state of the system state |^ n ) and the ancilla state |0), 

Win) |0) . (7) 

The quantum circuit then produces the following output state: 

X - (1 - ie^U) IVO |0> + \ (1 + ie*U) |VU |1) • (8) 

A projective measurement on the ancilla qubit in the computational basis {|0) , |1)} yields one of the two (unnormal- 
ized) states 

(l±ie*U)\il> in ), (9) 

which have mean energies either higher (for outcome |1)) or lower (for outcome |0)) than that of the initial state |^ n ). 
To justify this statement, let us expand the input state in the eigenvector basis {(e^)} of the Hamiltonian i7 s , 

\^in) =^2 c k |efc) . (10) 
k 

Note that the vector norms have the following form: 

\(l±ie^U)\e k )\ 2 = 2(l±sm0 k ) , (11) 

where (j>k = Ekt — 7 depends on the eigenenergy Ek of H s . In other words, each of the eigenvector weights \ck\ 2 is 
now multiplied by a factor 

l±sin0 fe (12) 

depending on the measurement result. To simplify our discussion, we will assume that one can always adjust the two 
parameters, 7 and t, such that 

- f < <\> k < f (13) 

for all fc's. Then, the factors (1 — sin^) are in descending order of the eigenenergies (see Fig. [1J), and the opposite 
is true for the factors (1 + sin^) . 

Therefore, apart from an overall normalization constant, the action of the operator (1 ± ie ll U) is to scale each of 
the probability weight \ck\ 2 by a factor of (1 ± sin (/>&); the probability weights scale to larger values, i.e., 

(l-sin^)/(l-sin^) > 1 (14) 

for the eigenenergy Ek < Ej in the cooling case (i.e., for outcome |0)), and vice versa for the heating case (i.e., for 
outcome |1)). 



Comments about the measurement 

The one-bit measurement result of the cooling module correlates with the energy of the state of the simulated 
system. After the measurement, the state of the system changes. If the state of the system is diagonal in the 
energy basis, this is strictly a Bayesian update. In general, this is nothing but a (partial) "collapse" of the wave 
function. If the measurement result is 0, the relative population of the low energy eigenstates of the system increases, 
Fig. [ljc). That is, if we post-select copies for which the results of a sequence of measurements have a sufficiently 
higher proportion of zeros than ones, then the average energy decreases, as in evaporative cooling. 

When each measurement corresponds to a short time evolution this technique is related to imaginary time evolution, 
as in diffusion Monte Carlo [? ]. If we do not post-select and the system is finite, the relative frequency of zeros and 
ones will peak on well defined values for different copies of the system. These well defined values are in a one-to-one 
correspondence with the possible energies of the Hamiltonian. The measurement of the energy can be viewed as a 
simulation of von Neumann's description of the quantum measurement process: the evolution of the system with 
Hamiltonian H couples with the position of a "measuring pointer." The pointer, in our case, is a binary ancilla qubit. 

Comments about "cooling" 

With our method, we do not result directly with a thermal state with Boltzmann distribution. Therefore, tem- 
perature is not well-defined. However, the entropy associated with the diagonal elements in the energy basis can be 
properly defined. Our method allows us to reduce the average energy and also the entropy. Of course, it is also 
possible to combine our method with other existing methods to project the exact thermal states, which will also allow 
us to define temperature. 

For example, if we start with the identity matrix / as our initial state, we can consider it as an infinite-temperature 
state. When we simulate a small time evolution and zero 7, then the scaling factor (e.g. see Fig. [T]) becomes 
approximately the same as the Boltzmann factor 

1 -sin<^ « e~ Ekt . (15) 

To further improve the fidelity, one may apply the projection method described in Ref. [? ? ]. We can repeat the 
same procedure to reach states with a lower temperature. 

Requirement for a fully scalable implementation 

In our experiment, we demonstrated a proof-of-principle experiment for implementing our cooling algorithm. How- 
ever, the experimental setup is specifically designed for cooling the polarization degrees of freedom. Here we discuss 
the requirements for a fully scalable implementation. First of all, we need a supply of fresh ancilla qubits prepared in 
the |0) state. It is not strictly required to have perfect |0) ancilla, but an error would decrease the cooling efficiency. 
Alternatively, if we are able to reset the ancilla qubit quickly, we can also employ only one ancilla qubit through 
out. Secondly, we will need to be able to maintain the coherence of the system qubits. Therefore, we expect error 
corrections will still be required, unless we can finish the cooling within the coherence time of the system. Thirdly, we 
will also need to have fast projective measurement on the ancilla qubits, which will allow us to be able to implement 
our feedback control system. An example is shown in Fig. [6] These requirements are not much different from those 
required for implementing a typical quantum algorithm. The main feature of this cooling algorithm is the modest 
requirement of using ancilla qubits, which could reduce the resources for expanding it into a fault-tolerant structure. 

2. SCALING ANALYSIS OF THE SIMULATED COOLING ALGORITHM 

Simplified model 

We study the complexity of preparing the ground state for easy of comparison with other methods. We focus on 
Hamiltonians with only two different energies. The generalization of the same problem to systems with more energies 
is straightforward. The degeneracy of the energy levels turns out not to be of any consequence for the resulting cost, 
only the projection probability on the lowest energy subspace is important. With that in mind, we will write the 
equations for a two level system to simplify the notation. 
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FIG. 6: An example of the scalable implementation of the quantum cooling algorithm. When the first step of cooling is imple- 
mented, a projective measurement (labeled by a question mark "?") is performed on the first ancilla qubit. The measurement 
result then determines if the second ancilla would interact with the first system or the second system. If cooling is resulted, 
then x = 1, and the same quantum circuit is applied to the same quantum system. It is same for both evaporative strategy and 
recycling strategy. However, if the measurement outcome is heating, then for the evaporative strategy, the system is rejected, 
and a new member is employed. For the recycling strategy, the same system is reset and sent to the same quantum circuit. 



We denote the two energy eigenstates by |eo) and |ei) with energies Eq and Ei, and define 

\^n) = Vp\ e o) + vT 3 ph> • (16) 

Let us ignore for the moment the possibility of filtering high energy states after too many "heating" measurements. 
The insight gained with this assumption can be applied to the more sophisticated algorithms in the text. We can then 
also assume, without loss of generality, that all the measurements are done together at the end of the computation. 
Therefore, if there are a total of k steps, we get the state (prior to the final measurement) 



E /(J) (V^ 1 / 2 + «) fe - J (l/2 - ^\eo)\Xj) + - p)(l/2 + 6)*-i(l/2 - ^| ei )|^)) , 



where and are properly normalized states, 



a = ^ sin(£ £ - 7) , & = ^ sin^it - 7) , 



(17) 



and j denotes the numbers of O's in both \Xj) and The parameter 7 is a controllable energy bias, see Fig. [I] in 

the main text. 

The probability to get j zeros on the final measurements of the ancillae is 

(1/2 + a) k -i(l/2 - aY + (1 - p) Q (1/2 + 6)^(1/2 - &)' . 

This is a mixture (with mixing probability p) of two binomial distributions. The left side binomial (corresponding to 
the ground state) concentrates on the interval 



and the right side binomial concentrates on the interval 



(18) 



(19) 



The first binomial corresponds to the probability of finding the ground state in the post-condition state after getting 
j zeros. Because we want to prepare the ground state, the binomials should have little intersection. This happens 
when 
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This is equivalent to 



Vk(b - a) > l - - (a 2 + b 2 ) . 



(20) 



The most relevant case is when the energies and the bias are small. We can do an expansion in terms of the small 
parameter At, where A is the gap and t is the evolution time. The parameters a and b are then small, O (At), using 
big-0 notation. We get 



\-{a 2 + b 2 ) = \+0{{At) 2 ) 



b — a 



At 



0{{At) 2 ). 



(21) 
(22) 



This gives the expected result 



At 



(23) 



That is, the number of measurements necessary to prepare the ground states has an scaling 



keO((At)~ 2 ) 



(24) 



If the binomials are separated, when we measure we get a number of zeros 

j ~ fc(l/2 - 6) , (25) 

which happens with probability 1 — p, and then we have collapsed on |ei), or we get a number of 

j~fc(l/2-o) , (26) 

which happens with probability p, and then we have collapsed on |eo). The number of measurements necessary to 
prepare the ground state scales like 

-^—2 ■ (27) 



p(At) 

Monte Carlo simulations 

We performed Monte Carlo simulations to test the scalings explained above for better algorithms, such as those 
implemented in the experiment. Notice that in the previous subsection we obtained scalings for the case when all 
the measurements are performed at the end or, equivalently, measurements are taken at each step but the result of 
the measurement is not taken into account at subsequent steps. Nevertheless, in the experiment adaptive action was 
taken taking into account the results of previous measurements (absorption and reflection boundaries, explained in 
main text). The intuition is that even an optimal adaptive action does not change the fundamental scaling. This is 
what we saw with Monte Carlo simulations. 

The actual experiment implemented an algorithm that keeps track of the difference x between the number of 
measurement results equal to and those equal to 1 in the ancilla, 

x = #0-#l . (28) 

The algorithm implemented in the experiment obeys two boundary conditions: 

1. If the number of measurement results equal to is less than the number of l's (x < 0), we reinitialize with the 
initial state. This is a "reflective" boundary condition. 

2. For some predefined constant c a b s , if x = c a b s then stop. This is an "absorbing" boundary condition. 

In addition, we need to define a constant Cbound such that we stop after reaching this number of steps. This is because 
this simple method can collapse to the excited state and never reach the absorbing wall (see previous subsection). 

We can analyze this algorithm following the intuition gained from the non-adaptive case explained above. Because 
we are interested in the asymptotic case with small energy differences, we set the energy bias 

7 = . (29) 

The position x of the random walker without boundary conditions was found to be distributed as a mixture of two 
(displaced) binomials distributions, centered around 

sm(Et) * k (30) 

with E taking the values of the ground and excited state energies, k the total number of measurements, and t is the 
evolution time. 

For the constant Cbound we can use our prediction on the number of steps needed to cool optimally 

Ci (awI ' (31) 

with a constant c\ that depends on the details of the algorithm. We want to show that this scaling holds for significant 
changes in At and p, when this parameters have small values. Because our process is not optimal, we actually use 

Cbound = 3 * pred (32) 

pred = d^ , (33) 
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FIG. 7: (a) Final fidelity for the ground state for the values explained in the text, (b) Final fidelity for the ground states for 
the values explained in the text when trowing away states that reach Cbound measurements, (c) Fidelity changes if we choose 
the reflective and absorbing boundary conditions according to the equations in the text, but using assumed gap values that 
do not corresponding to the real gap. The gap was kept at 0.02, but the gap used for the boundary conditions was the one 
indicated in the x-axis. 

where we use the constant c\ ~ 2.7 from the numerics of the optimal filtering algorithm (see below). Although in 
principle one can only use estimates of the gap A and the initial population p, we assume that this estimates are 
indeed good, and use the exact values in our numerics. Furthermore, we set the ground state energy Eq = — A and 
the excited state energy E\ = 0. These choices do not affect the fundamental scaling. 

Finally, we choose a value for the absorption boundary condition. We want this to be a value for x where states 
already have a good overlap with the ground state. We simply use the predicted average after the predicted number 
of measurements necessary to collapse 

Cabs = pred * s'm(E t) • (34) 

We performed Monte-Carlo numerics for p = 0.2 and At between 0.01 and 0.02, Fig.[7[ For these values, the number 
of maximum measurements (cbound) goes between 400000 and 100000, and the positions of the absorbing boundary 
condition goes between 1350 and 675. Nevertheless, the final fidelity with the ground state is quite constant, a clear 



indication that we understand the scaling correctly. The fidelity is always close to 0.76. For the error bars we use the 
expression 



A ±1.96^, (35) 

where a is the sample standard deviation, jl is the sample average, and we divide by the square root of the number 
of samples, V10000. 

As a simple improvement, we can remove from the calculation of the fidelity the states that reach the maximum 
number of steps, Cb OU nd- This corresponds to simply restarting again with an initial state. With the values chosen, 
only 20% of the measurement histories arrive to Cbound, so the probability of cooling with this method is also high. 
Further, the probability to always arrive at Cb OU nd decreases exponentially in the number of trials. The corresponding 
fidelities are depicted in Fig. [7| 

We also check that the fidelity does indeed vary substantially if we keep the gap fixed but change the reflective and 
absorbing bound conditions. This is plotted in Fig. [7] For these numerics the gap was kept constant at 0.02, but 
the reflective and absorbing boundary conditions were chosen as if the gap was different (between 0.01 and 0.02, as 
before). 

We also tested that the same fundamental scaling holds for an optimal refreshing schedule. The optimal refreshing 
schedule for the two level system is this: replace the current state with the initial state whenever the temperature of 
the current state is higher than the temperature of the initial state. This can be calculated with linear cost in the 
number of steps assuming that the energies of the ground state and the excited states are known, as well as the initial 
population of the ground state. In this case the boundary conditions and the bound on the maximum number of steps 



are unnecessary. Nevertheless, the distribution of the number of steps was found to be given by Eq. (31). The value 
of the constant c\ was found by fitting the Monte Carlo sampling data. 



Summary 



In summary, we focus on the complexity of preparing the ground state. The scaling is expressed in terms of the 
energy gap A between the ground state and the first excited state, the evolution time t of the simulated Hamiltonian 
for the weak energy measurement, and the population (probability) p of the ground state. The basic scaling is inversely 
proportional to (A t) 2 p. If possible, the theoretical optimal scaling occurs when t ~ 1/A, giving single-ancilla quantum 
phase estimation. Longer evolution times for each measurement require longer coherent times for the ancillae. It is 
also worth noting that the post-conditioned distribution is different for different evolution times and energy biases, 
which might be advantageous in some situations. Finally, in theory, the scaling with p can be improved to y/p using 
amplitude amplification . This comes at the expense of many more quantum operations, and therefore is less practical 
than the approach proposed here. 



3. FROM QUANTUM CIRCUIT TO EXPERIMENTAL IMPLEMENTATION 



In this section, we describe the experimental implementation of the quantum circuit for the cooling module in Fig.[T] 
We will focus on a single cooling module, and the cooling with respect to the Hamiltonian cr z . The <j x Hamiltonian 
is achieved by a rotation of the basis. Here we reproduce the relevant parts of the figures in the main text in Fig. [8] 
The goal here is to explain how to achieve a unitary transformation that is equivalent to the one in the quantum 
circuit shown in Fig. [8^. Our strategy is not to implement each operation separately, but to construct an effective 
transformation that produces the same overall result as the quantum circuit. 



Specification of certain simulation details 



1. In our experiment, for the Hamiltonian <j z , the mapping of the excited |e) and ground state \g) is |e) = \H) 
and \g) = |V), where \H) and \V) represent the horizontal and vertical polarization states of a single photon, 
respectively. For the a x Hamiltonian the excited state is \e) = (\H) + \V))/y/2 and the ground state is \g) = 
(\H)-\V))/V2. 

2. For the <j z Hamiltonian, the polarization state of the photon is mapped to the logical state of the system as 
|0) = \H) and |1) = \V). The two different path of the interferometer serves as the logical qubit of the ancilla. 
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FIG. 8: (Color online). The conversion from the quantum circuit to the experimental implementation. The abbreviations 
refers to polarization beam splitter (PBS), half- wave plate (HWP), quarter- wave plate (QWP), interference filter (IF), quartz 
plate compensator (PC), single photon detector (SPD), and fiber connector (FC). See the text for the explanation of the 
transformations . 



The system is prepared in the state a\0) + (3\1) (blue basis, see Fig[3]i) and the ancilla qubit is prepared as |0) 
(red basis) (step 1). 



3. In order to simulate cooling with the Hamiltonian cr x: the a x transform was applied to the input and output 
ports: we rotate the basis from \H) and \ V) to l/y/2(\H) + |V)) and l/y/2(\H) — \V)) This transform consists of 
a half- wave plate (HWP) with the angle setting at 22.5° and a tiltable quarter- wave plate (QWP). The logical 
circuit corresponding to the experimental setup is shown in Fig. [3^c), and the corresponding sequence of states 
is shown in Fig. [3^d). 

4. The input photon pulses for the simulation were generated by a laser with a central wavelength of 800 nm and 
a repetition rate of 76 MHz, attenuated to the single photon level. 



5. The measurement bases for the tomography were set by a QWP, a HWP and a polarization beam splitter (PBS). 
The photon was finally detected by a single photon detector (SPD) equipped with a 3 nm interference filter 
(IF). The probability of each output can be directly deduced from the single photon counts of the tomography 
measurements. 



Step-by-step instruction of the implementation of the cooling module 



1. In our setup, a photon is generated from the single photon source shown in Fig. [8J2. 

2. The photon is then passed through the state preparation module consisting of a polarization beam splitter 
(PBS), half- wave plate (HWP) and a quarter- wave plate (QWP). The resulting state is denoted as a\0) + /3|1) 
(blue basis, see Fig[8]i). The state associated with the path degrees of freedom is denoted as |0). Combining 
the two, we have the state 

(a|0)+/3|l})|0) (36) 
as the input state for the cooling module. This corresponds to step 1 of Fig. [8] 

3. The incident photon is then split by a polarization beam splitter, the horizontal polarized photons propagate 
directly while the vertical polarized photons reflect. This can be considered as a controlled-NOT operation. The 
resulting state becomes 

a|0)|0)+/?|l)|l) . (37) 

This corresponds to step 2 of Fig. [8j 

4. For each path, a unitary operation is realized by using two half- wave plates (HWP), and a quartz plate com- 
pensator (PC). The resulting state is 

a (a |0) + b |1)) |0) + (a |0) + b |1)) |1) . (38) 

This corresponds to step 3 and step 4 in Fig. [8] 

5. The two paths are then recombined on the same polarization beam splitter (PBS). The same controlled-NOT 
is applied, and the resulting is 

(aa |0) + 0b |1» |0) + (ab |1) + f3a |0» |1> . (39) 
This corresponds to step 5 in Fig. [8] 

6. For the path |1), a half- wave plate (controlled-not gate) is applied to the polarization state, giving the final state 

(aa\0) + 0&|1»|O) + (ab\0) + 0a|l))|l) . (40) 

Note that the unitaries are U \0) = U^l) = a\0) + where a = (1 + ie ie )/2 and b = (1 - ie ie )/2. This 
corresponds to step 6 in Fig. [8j 

7. At one output port, denoted as 0, the photon state is cooled and was sent to the next cooling module. At the 
other output port, denoted by 1, the photon is heated and was either rejected or recycled (re-prepared to the 
initial state and sent to the next cooling module). 



4. ADDITIONAL EXPERIMENTAL RESULTS 

Figure [9] shows the experimental results for the state transformation carried out by a single cooling module which 
simulates the Hamiltonian 

H s =a z . (41) 

The output state (photon) at the output port 1 (which corresponds to heating) is rejected. The ratio between \e) and 
\g) of the output state decreases by a factor 



(1 -sin0)/(l + sin0) 



(42) 



As a result, the output state becomes closer to the ground state with increasing 0, which is verified by our experimental 
results. For an energy bias angle 6 = 0° the state simply rotates in the equatorial plane of the Bloch sphere. For the 
input state 

(\e) + \g))/V2, (43) 

the output state as a function of the energy bias angle has Bloch sphere angles 

(a y ) = cos 6 and (a z ) = — sin 6 . (44) 

Figure |9|i shows the mean energies and ground state probabilities for different initial states (organized in panels 
with different colors) and different energy bias angles. The fidelities of the reconstructed states are shown in Fig. ^p. 
Red squares represent the experimental results and black dots represent the revised fidelities (see Methods Summary). 
The fidelities of the final states are larger than 0.978. Due to the compensation of symmetry noise, the revised fidelities 
are larger. 

Figure 10 shows the experimental results with evaporative and recycling strategies with initial input states 

(l/V5) |e) + (2/V5) \g) and (2/VS) \e) + (l/VS) \g) (45) 

respectively for the left and right panel, for Hamiltonian <j z (see also Fig. Qf-m)). Figure [lO^a) depicts the success 
probabilities at each step with evaporative and recycling strategies for input state (l/y/b)\e) + {2/y/b)\g) (0 — 10°). 
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FIG. 9: (Color online). Experimental results for a single cooling module with initial state (|e) + \g))/y/2 as a function of 
the energy bias angle 0. (a) Mean energies and ground-state probabilities of the reconstructed states for different energy bias 
angles. Different color panes represent the cases with different initial input states ("i" = a\e) + /3\g)) with relative populations 
\a\ 2 : |/3| 2 = 4 : 1, 1 : 1, 1 : 4. (b) The fidelities of the reconstructed states. Red squares represent the experimental fidelities. 
Black dots represent revised fidelities. Error bars are the corresponding standard deviations, (c) Tomography data. 
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FIG. 10: (Color online). Comparison of simulated cooling with evaporative and recycling strategies for the input states 
(l/\/5)|e) + (2/y/E)\g) (left) and (2/y/E)\e) + (l/y/E)\g) (right) for Hamiltonian a z . (a) Experimental results for the success 
probabilities in each cooling step with evaporative and recycling strategies for input state (l/y/b)\e) + (2/y/b)\g) (and = 10°). 
(b) Experimental results for the corresponding mean energies, (c-e) theoretical results for the mean energies for the evaporative 
strategy, the recycling strategy, and their difference, (f-h) Corresponding experimental results. Figures (i-p) show the same 
data but now with input state (2/y/E)\e) + (l/y/E)\g). 



Figure 10 ^b) shows the mean energies with evaporative and recycling strategies. Figures 10 [c-e) indicate the theoretical 
predictions for the mean energies for the cases with evaporative, with recycling, and their difference. Figures [lO^f-h) 



are the corresponding experimental results. The difference between the cases with evaporative and recycling is small 
compared to the case of Fig. [Z] Figure 101-p) shows the exact same plots for the input state (2/y/b)\e) + (\/y/h)\g). 



The difference between the cases with evaporative and recycling is larger than that in Fig. |10[c-h), as shown in 
Fig.^k-p). 

In the experiment, we performed three steps of cooling, we show the theoretical calculation for up to ten cooling 
steps in Fig. 11 The mean energies with evaporative are always smaller than the case with recycling. 




FIG. 11: (Color online). Theoretical results for the mean energies with ten cooling steps and Hamiltonian o z . Black squares 
represent the results with recycling and the red dots represent the case with evaporative. The initial state is (|e) + \g))/y/2 and 
the energy bias angle is set to 10°. 



